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The low-energy theory of graphene exhibits spontaneous chiral symmetry breaking due to pairing 
of quasiparticles and holes, corresponding to a semimetal-insulator transition at strong Coulomb 
coupling. We report a Lattice Monte Carlo study of the critical exponents of this transition as a 
function of the number of Dirac flavors N f , finding 5 = 1.25 ± 0.05 for N f = 0, 5 = 2.26 ± 0.06 
for Nf — 2 and 5 = 2.62 ± 0.11 for Nt = 4, with 7 ~ 1 throughout. We compare our results with 
recent analytical work for graphene and closely related systems, and discuss scenarios for the fate 
of the chiral transition at finite temperature and carrier density, an issue of relevance for upcoming 
experiments with suspended graphene samples. 
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Graphene, a single layer of carbon atoms ar- 
ranged in a honeycomb lattice [J provides a 
building block for more complex allotropes such as 
graphite (graphene sheets attached by van der Waals 
forces), fullcrenes (graphene spheres with pentago- 
nal dislocations) and nanotubes (cylindrically rolled-up 
graphene). In the absence of electron-electron interac- 
tions, the valence and conduction bands of graphene are 
connected by two inequivalent "Dirac points" , around 
which the low-energy excitations are massless quasiparti- 
cles with a linear dispersion relation and a Fermi velocity 
of v F ~ c/300 Hi]. 

Such a scmimctallic band structure 
is, unfortunately, unsuitable for many electronic applica- 
tions, which depend crucially on the ability to externally 
modify the conduction properties, as routinely done with 
semiconducting devices. The quest to engineer a bandgap 
in graphene has thus been propelled to the forefront of 
current research. Hitherto suggested solutions include 
gap formation due to interaction with a substrate [1] , in- 
duction of strain [f| , and geometric confinement by means 
of nanoribbons or quantum dots [6j. 

The low v F in graphene indicates that the analog of 
the fine-structure constant of Quantum Electrodynam- 
ics (QED) is « 9 ~ 1, and thus the Coulomb attrac- 
tion between electrons and holes may play a significant 
role in defining the ground-state properties. An intrigu- 
ing possibility is that spontaneous formation of excitons 
(electron-hole bound states) and the concomitant break- 
ing of chiral symmetry may turn graphene into a Mott 
insulator. While the strength of the Coulomb interac- 
tion precludes a perturbative approach, previous (ap- 
proximate) analytic studies [7| at the neutral point (zero 
carrier density n) and zero temperature T have addressed 
the appearance of an excitonic gap as a function of a g 
(see Fig. [T]). Such treatments suggest that the transition 
into the insulating phase should be governed by essential 
singularities rather than power laws, a behavior known 
as Miransky scaling 

In our recent Lattice Monte Carlo (LMC) study @ , in- 
dications were found that the chiral transition is of second 



order, with well-defined critical exponents. Subsequently, 
in Ref. [l(| we provided a rough estimate of the critical 
exponents as S ~ 2.3, j3 m ~ 0.8 and 7 ~ 1, although 
a more precise determination was not possible due to 
insufficient data on large enough lattices. Nevertheless, 
Miransky scaling and classical mean-field exponents were 
found to be disfavored. 

The aim of the present work is to provide a more rig- 
orous and comprehensive determination of the quantum 
critical properties for N * = 0, 2 and 4 Dirac flavors, as 
well as to contrast these results with recent analytical 
and simulational work for graphene and related theories. 
We also briefly elaborate on the mechanisms that inhibit 
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FIG. 1: (Color online) Phase diagram in the (Nf and /3) 
plane. The gapped phase is bounded by a critical Coulomb 
coupling a c = (47r/3 c )~ 1 and a critical number of fermion fla- 
vors Nf c . The coupling on a Si02 substrate is denoted by 
a Si0 2 ' an d f° r suspended graphene by a susp . Inset: hypo- 
thetical phase diagram in the (n, T) plane. At low T, sus- 
pended graphene exhibits semimetallic properties whenever 
the carrier density n exceeds a characteristic value n* . At 
the neutral point the semiconducting behavior persists un- 
til T ~ ejp = hv F \/-Kn* , where the transition may be of 
Berezinskii-Kosterlitz-Thouless type or a crossover. 
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exciton formation at non-zero T and n, and their con- 
nection to other systems. 

The LMC studies of Refs. [I OH ED] suggest that the 
low-energy theory of graphene is an appropriate starting 
point for a quantitative analysis. This is defined by the 
Euclidean action 



/ 



d 2 xdt^ a D[A ]^ a 



2e 2 



d 3 xdt (d % A ) 



(1) 



with the Dirac operator 

D[A ] = 7o (d + iA ) + v 7i d i + m t, i = l,2 (2) 

where the ip a with a = 1 , . . . , N j are four-component 
spinors in 2+1 dimensions, A is a Coulomb field in 3+1 
dimensions, and the case of a graphene monolayer is re- 
covered for Nj — 2 in the limit m — + 0. Furthermore, 
a g = e 2 /(AitvEq) with the inverse coupling (3 = ve /e , 
such that screening by a substrate is reflected in the di- 
electric constant e . 

The gauge term of E q fl] ) is discretized in the non- 
compact formulation lldj|T The staggered discretiza- 
tion [l2j of the fermionic component of Eq. fl} is pre- 
ferred, as chiral symmetry is then partially retained at 
finite lattice spacing. As N staggered flavors correspond 
to N f — 2N continuum Dirac flavors UM, the case of 



V 

= 2 is recovered for N — 1, 



giving 



Se[x,X,U ] = 
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where the x n are staggered fermion spinors, and the site 
indices (m, n) are restricted to a 2+1 dimensional sublat- 
tice. Invariance under spatially uniform, time-dependent 
gauge transformations is retained by the link variables 



U, 



= u„ 



exp(j# n ), where 9 n is the lattice gauge 



field. For v = 1, the staggered form of Eq. ([2]) is 



(4) 
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where n = (— l) n ° and r] 2n — (— l) n o+™i. Our simula- 
tions use the Hybrid Monte Carlo (HMC) algorithm with 
N pseudofermion flavors on a 2+1 dimensional space- 
time lattice of extent L, such that 9 also propagates in 
the third spatial dimension of extent L z . Further details 
are given in Refs. [l(J 0] • 

We now seek to characterize the critical exponents 
of the chiral transition in graphene. The spontaneous 
breakdown of chiral symmetry in Eq. ([1]) is signaled by 
a non-zero condensate a = (xx)- The mass term in 
Eq. © breaks chiral symmetry explicitly, generating a 
non-vanishing condensate, which is otherwise not possi- 
ble at finite volume. The appearance of a gap in the 
quasiparticle spectrum of graphene at a critical coupling 
[3 C is then marked by a ^ for m — > 0. However, the 



"chiral limit" m — > cannot be approached directly, as 
that limit corresponds to a very large fermionic corre- 
lation length, especially in the vicinity of [3 C due to the 
appearance of Goldstone modes. Practical simulations 
are performed at finite m , such that the limit m —> 
is reached by extrapolation, for which it is useful to also 
study the susceptibility Xi = dcr/dm and the logarith- 
mic derivative R = dhia/dhim . An instructive way 
to determine /3 C and the critical exponents is by fitting 
an equation of state (EOS) m = f(a,j3) to simulation 
data at finite m . Knowledge of f(cr,/3) with good pre- 
cision close to the transition then allows for an educated 
extrapolation to the chiral limit. 

We have considered the EOS successfully applied [l5| 
to lattice QED, 



m X(/3) = rMAW+^a) 



(5) 



where X(J3) and Y{j3) are expanded around j3 c such that 
X(p) = Xo+X^l-P/fiJ andr(/3) = Y^l-f}/^). The 
dependence on a is given by = a h and / 3 (er) = a s , 

where b = 5 — l//3 m . The critical exponents are 
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FIG. 2: (Color online) Chiral condensate a (left panel) and 
susceptibility Xi (right panel) for Nf — 2. Data for L = 
28, L z = 8 are indicated by squares, and for L — 32, L z = 
12 by circles. The lines represent a \ 2 fit to rr and and 
extrapolation m — > using Eq. lj5}. The open datapoints are 
excluded due to finite-volume or lattice spacing effects. The 
optimal fit is (3 C = 0.0738 ± 0.0010 and 5 = 2.23 ± 0.06, with 
X = 0.36 ± 0.05, X 1 = -0.13 ± 0.02 and Y 1 = -0.15 ± 0.02. 
The errors are of statistical origin. The method of analysis is 
described in detail in Refs. [gt fioj] - 
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FIG. 3: (Color online) Chiral condensate a (left panel) and 
susceptibility \i (right panel) for Nt = 4, for lattices up to 
L = 24, L z = 8. Inset: logarithmic derivative R for different 
p' 1 . The optimal fit is /3 C = 0.0499 ± 0.0010 and S = 2.62 ± 
0.11, with X = 0.19 ± 0.05, X x = -0.09 ± 0.02 and Y x = 
-0.08 ± 0.02. See also Fig. H 



FIG. 4: (Color online) Chiral condensate a (left panel) and 
susceptibility \i (right panel) for Nf = 0, with L — 36, L z = 
20 (squares) and L — 48, L z = 20 (circles). The optimal fit, 
with fixed 5 = 1.25 is (i c = 0.158 ± 0.001, X = 0.16 ± 0.02, 
X 1 = -0.10 ± 0.05 and Y 1 = -0.11 ± 0.02. The range S = 
1.25 ± 0.05 yields /3 C = 0.16 ± 0.02. 



which are assumed to obey the hyperscaling relation 
(3 m (5 — 1) = 7. It should also be noted that R — > 1/5 
for to — > at (3 — (3 C . Our results and analyses in terms 
of Eq. ([5]) are shown in Fig. [2] for Nf = 2, in Fig. [3] for 
Nf = 4, and for the quenched case Nf = in Fig. 0J All 
of our results are consistent with b = 1.00 ± 0.05, hence 
we conclude that 7 ~ 1 based on the hyperscaling rela- 
tion, such that the remaining exponent to determine is 6. 
Based on the EOS analysis and the logarithmic deriva- 
tive R (see Fig. E]), we find 6 = 2.26 ± 0.06 for N f = 2, 
6 = 2.62 ± 0.11 for N f = 4, and 5 = 1.25 ± 0.05 for 
Nj = 0. We observe that finite volume effects decrease 
with increasing Nf and that datapoints for small (3 and 
large m in the broken phase are not well described by 
Eq. ([5]), likely due to a small correlation length associated 
with a growing excitonic gap. 

An increase in 8 with is consistent with the LMC 
results of Ref. [lH where a similar trend was found, cul- 
minating at 3.6 < S < 6 for Ay = N^ c ~ 4.8 where 
the chiral transition disappears. Such behavior is remi- 
niscent of the Thirring model in 2 + 1 dimensions 
where S ~ 2.8 for N* — 2, reaching S ~ 7 at a criti- 
cal flavor number of N i c ~ 6.6. Extensive LMC studies 
of QED have found S ~ 2.2 for N f = [13, while for 
QED with dynamical fermions (5 ~ 3 [l^|. The case of 
QED in 2 + 1 dimensions (QED 3 ) is noteworthy as the 
LMC study of Ref. [3 yielded 5 ~ 2.3 for N f = 1 and 
S ~ 2.7 for ^ = 4, which are suggestive of our val- 
ues for graphene, although spontaneous chiral symmetry 
breaking in QED 3 is difficult to establish as the order 
parameter can be exponentially suppressed for large Nf. 

The gap-equation analysis of Ref. [7| reported j3 c ~ 
0.16 for N f = and (3 C ~ 0.066 for N f = 2, which are 
in qualitative agreement with our results. However, the 



transition of Ref. [7| is of infinite order and vanishes for 
Nf = 4. These discrepancies are smallest for A^ = 0, 
where our results approach 6 = 1. Our observations 
are thus in line with indications [l9| that the critical 
exponents, as obtained from Schwinger-Dyson equation 
(SDE) analyses, may be dependent on the chosen resum- 
mation scheme. 

An effective theory containing both the order param- 
eter and the Dirac quasiparticles as dynamical fields has 
recently been developed in Ref. [2(| ■ Based on an expan- 
sion to leading order around e = 3 — d spatial dimensions, 
the long-range ~ 1/r Coulomb tail was found to be ir- 
relevant in the renormalization group (RG) sense, such 
that the chiral transition could then be described using 
only short-range interactions of the Gross-Neveu- Yukawa 
form, yielding the estimates 7 ~ 1.25 and S ~ 2.8 for the 
critical exponents at Nf = 2, in qualitative agreement 
with our present findings, as well as with large- ./V^ calcu- 
lations of the RG flow [2l|. However, our results are not 
compatible with S = 2 + 0(l/N f ) found in Refs. [USl, 
which is surprising as the Gross-Neveu- Yukawa theory is 
expected to interpolate between S = 19/5 ~ 4 at Nf = 
and 5 = 2 in the Nf — » 00 limit [20j . It is not clear, as no 
chiral transition exists in the graphene theory above the 
critical flavor number Nf c = 4.8 [ll(, how to consistently 
compare our results with large-Aj estimates. 

What is the fate of the semimetal-insulator transition 
at non-zero temperature? On the basis of the Mermin- 
Wagner theorem [23| , one expects either a crossover or a 
Berezinskii-Kosterlitz-Thouless (BKT) transition [3] at 
a critical temperature T c . The most compelling exper- 
imental evidence so far for a BKT transition has been 
reported in Ref. [25[, where graphene samples on a sub- 
strate were subjected to transverse magnetic fields up to 
B ~ 30 T. In the temperature range of 10 K to 1 K, a 
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growth in the resistivity by a factor of ~ 200 was ob- 
served, and attributed to the "magnetic catalysis" pre- 
dicted in Refs. [H, H3|. At B = 0, the resistivity of 
annealed suspended graphene was observed (28j to in- 
crease by a factor of ~ 3 over a temperature range of 
200 K to 50 K, while also changing character from metal- 
lic to semiconducting. A study of the low-energy theory 
of graphene at non-zero T is thus clearly called for, pos- 
sibly along the lines of Ref. [2i|, which considered the 
Gross-Neveu model in 2 + 1 dimensions. 

At low T, the large extent of the imaginary time di- 
mension renders the system effectively three-dimensional, 
such that a chiral transition may still be observed at a 
critical density n*. Interestingly, away from the neutral 
(unpolarizcd) point, non-relativistic Fermi systems (such 
as the asymmetric Fermi liquid in the context of ultracold 




) can undergo tran- 
before reverting to 



atoms and dilute neutron matter 
sitions into exotic phases [3l|, [32 
a fully polarized normal state. Whether the low-energy 
theory of graphene exhibits such phenomena is currently 
unknown. 
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Nf = indicates that S ~ 1.2, while the case of Nf = 2 is 
compatible with 8 ~ 2.2. Finite-volume effects are substantial 
for Nf — 0, at large (3 and small m . 
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